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1.  INTRODUCTION 


Ballistic  experiments,  like  most  other  experiments,  usually  have  random  errors  associated 
with  measurements.  The  errors  come  from  many  factors,  such  as  differences  in  operational 
procedures,  environmental  conditions,  and  material  properties.  Of  these,  operational  procedures 
can  usually  be  maintained  constant.  However,  there  are  uncertainties  associated  with  the  measured 
variables,  such  as  yaw,  pitch,  impact  velocity,  etc.,  which  can  affect  the  response  variables,  e.g., 
residual  mass,  residual  velocity,  etc.  In  addition,  there  is  little  that  can  be  done  to  control 
environmental  conditions  except  to  monitor  the  weather  conditions. 

An  additional  set  of  measured  variaUes  that  are  generally  not  treated  well  in  ballistic 
experiments  is  the  material  properties.  In  most  cases,  absence  of  detailed  descriptions  about  the 
properties  adds  a  source  of  uncontrolled  error  in  the  response  variable.  In  some  cases,  material 
property  variation  could  be  the  primary  factor,  since,  in  general,  materials  are  neither  homogeneous 
nor  isotropic  but  vary  from  piece  to  piece  and  point  to  point  within  a  sample.  These  departures  from 
homogeneity  may  be  the  result  of  local  chemistry  variations,  microstructure  variation,  textures 
developed  during  processing,  or  other  factors.  Each  definitely  contributes  to  the  scattering  in  the 
experimental  results. 

The  purpose  of  this  report  is  to  investigate  the  way  variations  in  material  properties  might 
affect  the  outcome  of  a  series  of  ballistic  penetration  computations.  As  currently  configured, 
Eulerian  as  well  as  Lagrangian  computational  codes  suitable  for  penetration  computations  assume 
that  the  materials  are  homogeneous  and  isotropic.  Additionally,  other  than  stress/strain  criteria  for 
failure,  there  are  no  other  explicit  failure  criteria  built  in  simulation  packages  such  as  HULL.  A 
number  of  Government  laboratories  and  universities  are  striving  to  model  such  failure  mechanisms 
as  adiabatic  shear  failure  (Hauver  et  al.  1992).  These  models  start  at  the  most  fundamental 
microscopic  level  (for  instance,  the  inclusion  of  the  “shear  band”  failure  mechanism  in  the  code); 
their  incorporation  in  a  large  computational  code  that  can  run  to  completion  within  memory  and 
time  constraints  is  not  a  certainty. 
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This  study  attempts  to  attack  the  problem  from  the  opposite  end,  zeroth-order  macroscopic 
direction.  It  is  a  first  attempt  to  add  stochastic  variations  to  the  material  property  subroutines  in  an 
otherwise  deterministic  simulation  package.  The  hope  is  to  gain  insight  on  the  level  of  modeling 
required  for  accurate  prediction  of  ballistic  results  and  to  provide  researchers  with  a  direction  for 
further  refinement. 

As  a  sample  problem,  the  relatively  simple  case  of  a  long  rod  penetrator  striking  an  oblique, 
finite-thickness  plate  was  computed  with  and  without  stochastic  variations  in  material  properties. 
The  results  are  presented  and  discussed. 


2.  METHOD 


The  HULL  eode  handles  material  failure  by  giving  the  user  the  option  to  choose  from  the 
following  three  different  failure  criteria;  (1)  the  maximum  principal  stress  criterion  for  spall  failure 
in  plane  strain  (FAIL=1  in  the  HULL  input  deck);  (2^  the  maximum  principal  strain  criterion  added 
for  failure  in  plane  stress  or  ductile  failure  (FAIL=1  STRAIN=1)  (in  the  present  study,  both  ultimate 
stress  and  strain  were  used  to  initiate  material  failure);  (3)  the  triaxial  states,  based  on  a  PA'  model 
(a  material  failure  domain  to  strain-stress  relationship),  used  to  determine  whether  fracture  can 
occur.  The  third  is  most  complex  and  debatably  more  realistic  in  behavior  but  requires  extensive 
material  property  data  (FAIL=2  STRAIN=1)  (Matuska,  Osborn,  and  Pibum  1991).  All  these  failure 
models  result  in  material  separation  or  void  inclusion.  In  HULL,  the  most  numerically  reliable 
method  of  simulating  these  conditions  has  been  to  “inject”  air  into  the  calculation  as  the  “void” 
material. 

During  the  initial  stage  of  this  study,  the  first  failure  model  was  used.  Since  this  model  was 
elementary,  the  second  model  was  eventually  employed.  The  third  model  uses  more  detailed 
description,  but  since  extensive  material  property  data  were  requited,  it  was  not  possible  to 
implement  at  this  time.  Consequently,  this  study  proceeded  primarily  with  the  second  failure 
model. 


2 


Careful  study  of  the  simulation  package  HULL  revealed  that  the  material  failure  mechanism 
is  described  in  subroutines  “hydro”  (in  subprogram  huh)  and  “mgrun"  (in  sul^rogram  eos),  where 
the  stress  and  strain  are  hist  calculated  in  the  usual  deterministic  way  and  then  compared  with  the 
ultimate  stress  or  strain  to  determine  whether  the  particular  computational  material  ceU  fails.  The 
chaUenge  was  to  find  a  way  to  add  stochastic  behavior  into  this  otherwise  deterministic  failure 
mechanism. 

With  a  single  specimen  of  material  undergoing  a  uniaxial  strength  test,  there  is  the  familiar 
stress-strain  relationship  as  shown  in  Figure  1.  Here  the  material  stretches  linear-elasticaUy  until  it 
reaches  a  certain  yield  stress.  After  this  point,  the  relationship  is  no  longer  elastic.  The  material 
continues  to  elongate  until  it  reaches  a  point  where  it  can  no  longer  sustain  an  increase  in  stress,  and 
there  it  breaks  or  fails.  For  a  single  test,  one  obtains  a  smooth  stress-strain  curve.  But  the  materials 
are  not  leaUy  uniform;  they  are  not  the  same  from  sample  to  sample.  The  differences  may  be  caused 
by  chemical  inhomogenieties,  texturing  caused  by  mechanical  working  during  processing,  and  so 
forth.  Consequently,  there  will  be  variations  in  this  curve,  and  the  location  of  the  failure  point  is  not 
the  same  from  specimen  to  specimea 

To  quantify  this  variation,  a  variable  called  the  “probability  of  failure”  as  a  function  of  stress 
was  used.  This  implies  that  there  is  a  certain  probability  of  failure  associated  with  each  stress  level. 
The  material  does  not  fail  immediately  wlwn  a  certain  ultimate  stress  is  reached,  but  rather  the 
probability  of  failure  increases  around  that  stress  level  Oy  at  certain  cell  location  ij  within  an  interval 
6,  as  shown  in  Figure  2.  A  certain  statistical  distribution  of  this  failure  behavior  is  assumed,  where 
the  probability  of  failure  ranges  from  0  (i.e.,  no  failure)  to  1  (i.e.,  certain  failure).  Wth  this 
statistical  distribution,  it  is  possible  to  bring  some  stochastic  behavior  into  the  otherwise 
deterministic  failure  mechanism  with  the  following  steps: 

(1)  calculate  the  stress  in  a  particular  material  cell  being  looked  at,  sequentially  through  all 
cells  in  all  coordinates, 

(2)  determine  the  probability  of  failure  from  the  statistical  distribution  of  failure  as  a  function 
of  stress. 
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Figure  1.  The  ordinary  stress-strain  relationship. 


Figure  2.  Probability  of  failure  as  a  ftmction  of  stress. 
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(3)  draw  a  random  number  from  a  uniform  random  number  (URN)  routine  (the  number  falls 
uniformly  between  0  and  1), 

(4)  compare  the  probability  of  failure  with  the  random  number:  if  URN  <=  probability  of 
failure,  then  the  material  cell  fails. 

In  the  stochastic  treatment  of  the  relationship  between  stress  and  material  failure,  the  amount 
of  stress  necessary  to  cause  such  a  failure  becomes  a  random  variable  which  is  assumed  to  be 
normally  distributed  with  an  expected  value  equal  to  the  original  stress  threshold.  Other 
distributions  can  be  applied  if  there  are  logical  reasons  to  do  so;  however,  in  this  case  the  normal 
assumption  seems  intuitively  applicable.  The  nonnal  probability  density  function  is  tlM  familiar 
bell-shaped  curve;  its  cumulative  distributitm  function  results  from  the  integration  of  this 
probability  density  function  (Walpole  and  Myers  1978).  For  this  procedure,  the  cumulative 
distributitm  function  was  truncated  so  that  it  covered  a  sufficient  range  of  stress  (±5  in  Figure  2) 
without  reaching  infinite  values.  This  truncated  normal  distribution  was  then  used  as  the  functional 
relationship  between  the  amount  of  stress  and  the  probability  of  material  failure  in  all  the 
computational  cells. 

A  uniform  random  number  is  chosen  from  a  routine  which  uses  the  computer  system  time  as 
the  seed  for  initiation  of  calculation.  This  is  a  traditional  technique  used  in  Monte-Carlo 
simulations.  It  provides  equally  likely  numbers  between  0  and  1  which  can  be  compared  with  some 
cumulative  distribution  function.  The  computation  follows  the  previously  mentioned  steps  to 
decide  whether  the  material  in  a  particular  computational  cell  fails.  If  it  does,  then  a  void  is  inserted 
into  that  cell.  In  this  example,  if  the  probability  of  material  failure  from  the  cumulative  normal 
distribution  function  is  0.9,  then  approximately  9  times  out  of  10  the  uniform  random  number  will 
be  less  than  or  equal  to  the  probability  of  failure.  Thus,  in  10  replications  of  the  model,  there  will 
be  ^proximately  9  failures.  To  ensure  sufficient  stochastic  effects,  the  penetration  problem  was 
repeated  100  times  on  a  Cray-2  computer,  calculating  out  to  200  ps  physical  time.  For  each 
replication,  the  final  locations  of  both  the  tip  and  the  tail  of  the  penetrator  were  determined,  and  the 
resulting  lengths  were  statistically  analyzed. 
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An  additional  algorithm  to  add  stochastic  behavior  to  the  strain  fail  criterion  has  been 
formulated  and  is  similar  to  the  one  for  the  stress  fail  mentioned  previously. 

Later,  some  modifications  were  made.  This  wUl  be  described  in  the  next  section. 

3.  RESULTS 

The  problem  chosen  for  the  statistical  study  is  a  tungsten  alloy  penetrator  against  a  high 
hardness  annor  (HHA)  plate  (Magness,  Farrand,  Rensselaer  1986).  An  HHA  plate  of  3.175-cm 
(1.25  in)  thickness  at  60*  obliquity  was  tested  against  a  tungsten  alloy  penetrator  (radius=0.414  cm, 
length=10.26  cm,  L/I>=  12.39)  at  1.1  km/s  (see  Appendix  for  the  input  deck). 

The  initial  configuration  is  shown  in  Figure  3. 
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Figure  3.  Armor  example  at  0  ps. 


The  computation  created  data  dumps  eveiy  50  ps,  until  200  ps  was  reached.  These  were 
plotted  and  are  displayed  in  Figures  4-7. 
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Figure  4.  Annor  example  at  SO  tts. 
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Figure  S.  Aimor  example  at  100  ps. 
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Figure  6.  Annor  example  at  150  ps. 
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Figure  7.  Arnior  example  at  200  ps. 


For  each  of  the  three  percent  deviations  associated  with  the  cumulative  distribution,  namely 
ultimate  stress  ±S%,  ±10%,  and  ±15%  (±5  in  Figure  2),  the  calculation  was  carried  out  100  times. 
A  calculation  with  0%  deviation,  i.e.,  without  any  statistical  randomness,  was  also  carried  out  to 
serve  as  a  basis  for  comparison.  The  final  lengths  of  the  penetrator  after  200  ps  were  collected  and 
statistically  processed,  as  displayed  in  Table  1. 
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Table  1:  Final  Lengths  vs.  Initial  Randomness  (First  Calculation) 


Percent  Dev. 

of  the 
cumulative 
normal 
distribution 

Median 
of  the  100 
computed 
final  lengths 
(cm) 

btw  stations 

1  and  2 

Standard 

Dev. 

of  the  100 
computed 
final  lengths 

0.00 

2.2488 

0.05 

2.24560 

0.074155 

2.23710 

0.078487 

2.24365 

0.073094 

0.10 

2.28730 

0.076906 

2.27935 

0.078190 

2.28375 

0.074500 

0.15 

2.29045 

0.071916 

2.28455 

0.077485 

2.29350 

0.079646 

Notice  that,  for  each  percent  deviation  of  the  cumulative  normal  distribution,  the  patched 
HULL  has  been  run  100  times.  But,  in  order  to  see  also  the  statistical  character  of  this  computation, 
this  has  been  done  three  times  for  each  of  the  percent  deviations,  as  displayed  in  the  three  rows  for 
each  percent  deviation  in  Table  1 .  The  results  indicate  a  distribution  of  final  penetrator  lengths,  and 
this  randomness  shows  a  mild  trend  corresponding  to  the  trend  in  the  percent  deviation  of  the 
cumulative  normal  distribution  between  the  probability  of  failure  and  stress/strain.  The  in-group 
randomness  (last  column  in  Table  1)  remains  about  constant. 

After  some  deliberation  on  these  results,  it  was  noticed  that,  while  the  seed  (which  is  required 
to  initiate  the  calculation  for  a  random  number)  used  in  the  previous  computation  did  change 
randomly,  which  is  desirable,  the  randomness  was  more  than  necessary.  While  the  computation 
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nins  sequentially  in  the  spacial  coordinates,  when  the  computation  sweeps  into  the  next  time  step, 
the  seed  used  in  each  material  cell  may  not  be  exactly  the  same  seed  used  during  the  previous  time 
step  at  exactly  that  same  location.  Tliis  variation  along  the  time  axis  did  not  seem  to  correctly 
describe  the  variation  in  the  material  property.  In  effect,  this  numerical  process  tended  to 
homogenize  the  material  properties.  Therefore,  the  next  attempt  removed  this  time-dependent 
variation. 

One  way  to  implement  the  randomness  into  the  computation  without  a  time-dependent 
randomness  is  to  reset  the  seed  every  time  step  at  the  same  material  cell  to  be  exactly  the  same  as 
the  value  at  the  beginning.  The  computation  results  from  this  formulation  are  disfdayed  in  Tible  2. 


Table  2:  Final  Lengths  vs.  Initial  Randonmess  (Second  Calculation) 


Percent  Dev. 

of  the 
cumulative 
normal 
distribution 

Median 
of  the  100 
computed 
final  lengths 
(cm) 

btw  stations 

1  and  2 

Standard 

Dev. 

of  the  100 
computed 
final  lengths 

0.00 

4.1682 

0.05 

4.18219 

0.02968 

0.10 

4.18715 

0.01083 

0.15 

4.20477 

0.00793 

(Note;  The  grid  and  station  locations  have  been  modified  from  Tbble  1  to  Table  2  during  the 
evolution  of  the  study,  such  that  the  final  lengths  between  stations  show  different  values.  But  the 
statistical  meaning  of  the  study  is  not  changed.) 


The  results  in  Table  2  show  some  trend  relationship  between  the  randomness  in  the  final 
penetrator  length  and  the  randomness  inherent  in  cumulative  normal  distribution  fimction.  Since 


13 


higher  percent  deviation  is  likely  associated  with  higher  fail  rate  of  the  target  material,  thus,  at  lower 
percent  deviations,  there  is  less  penetration  into  the  target  and  more  compression  of  the  penetrator, 
resulting  in  a  smaller  final  length.  So  the  relationship  between  the  randomness  in  the  failure  model 
and  the  randomness  in  the  final  length  seems  to  be  a  result  of  the  compression  of  the  penetrator 
rather  than  a  result  of  its  erosion. 

Furthermore,  since  the  difference  is  so  small,  the  randomness  itself  may  be  covered  by  other 
noise  such  as  that  from  the  computational  accuracy.  In  addition,  the  in-group  randomness  (last 
column)  shows  a  clearer  trend  than  that  in  Ihble  1. 

The  differences  in  the  median  final  lengths  in  Table  2  from  those  in  Table  1  reflect  the  changes 
that  occurred  in  the  input  geometry  (grid  arrangement  and  station  locations)  during  the  evolution  of 
development  (see  Appendix).  These  difference  are  not  significant  toward  the  understanding  of  the 
statistical  behavior  of  the  material  strength. 


4.  DISCUSSIONS 

Two  different  methods  to  add  stochastic  behavior  to  the  deterministic  simulation  package 
HULL  were  investigated. 

The  results  from  the  first  method,  continued  regeneration  of  random  numbers  with  a  new  seed 
at  each  time  step,  showed  certain  relationship  between  the  randomness  associated  with  the 
cumulative  normal  distribution  function  and  the  randomness  in  the  final  length  of  the  penetrator. 
However,  the  varying  seed  in  each  cell  in  time  is  more  than  necessary.  A  reasonable  fix  was  to 
remove  the  variation  in  time. 

The  results  from  the  second  method,  resetting  the  seed  value  to  its  initial  value  for  all 
subsequent  time  steps,  showed  a  mild  relationship,  which  seemed  to  be  driven  mote  by  the 
compression  of  the  penetrator  than  by  its  erosion.  In  addition,  there  is  also  a  trend  in  the  in-group 
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randomness  Oast  column)  relative  to  the  randomness  in  the  failure  model.  This  trend  was  small  and 
possibly  tainted  by  faint  computational  noises. 


Further  study  in  the  last  approach  has  revealed  the  fact  that  the  computational  grid  in  an 
Eulerian  formulation  does  not  follow  the  material.  In  other  words,  the  Eulerian  grid  does  not 
capture  the  material  property  in  the  same  cell.  In  an  Eulerian  code  such  as  HULL,  the  calculation 
sweeps  through  grid  nodes  in  the  Eulerian  space.  Modifications  to  the  code  change  the  calculation 
at  each  grid  point  in  the  Eulerian  grid  but  not  the  material.  So  this  formulation  does  not  provide  a 
random  material  property  at  the  same  material  location,  which  is  unfortunately  essential  for  the 
study  of  the  material  failure. 

In  order  to  modify  the  calculation  in  association  to  each  certain  material  particle,  a  more 
sophisticated  way  must  be  formulated,  such  that  the  calculation  will  follow  the  material  instead  of 
merely  the  Eulerian  grid.  Therefore,  the  next  reasonable  approach  is  to  find  a  way  to  add  the 
stochastic  variation  in  a  Lagrangian  fashion,  which  is  not  quite  as  straightforward  in  an  Eulerian 
code.  Furthermore,  Lagrangian  codes  usually  have  limitations  in  dealing  with  failure  problems. 


However,  since  Table  1  (first  method,  which  has  more  randomness  than  necessary)  and 
Table  2  (second  method)  both  showed  a  trend  between  the  final  lengths  and  the  initial  randomnesses 
in  the  model,  it  is  likely  a  similar  trend  will  also  be  present  in  a  better  formulated  calculation  based 
upon  a  Lagrangian  grid. 

This  study  exposed  some  inadequacies  in  two  different  approaches  and  provided  iUumination 
for  future  study  in  the  statistical  analysis  of  material  strength. 
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APPENDIX: 

TYPICAL  INPUT  DECKS 
FOR  THE  FIRST  AND  SECOND  CALCULATIONS 
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A- 1 .  INPUT  DECK  FOR  THE  FIRST  CALCULATION 


A  typical  iiqnit  file  to  initiate  the  HULL  nin,  called  keel.in,  follows: 


keel  prob  00221.0100 
diinen«2  geom-1 
imax-50  jmax-100 
fail=l  strain=l 
nm=3  air“l  hha“2  walloy-3 
visc=l  d.visc-2 
nop=2  nstn-2 
header 

walloy  -  hha  @  1.1  km/s,  stochastic  stress+strain  failure 

mesh 

x0=-5  xmax=5  y0=-ll  ymax=9 
generate 

package  walloy  v»l.le5 

rectangle  xl=-.414  x2=.414  yl=-10.26  y2=-.414 
circle  xc»0 .  yc=-0.414  radius=0.414 
package  walloy  v-l.leS 

circle  xc=0 .  yc=-0.414  radius=0.414 
package  hha 

rectangle  xl=-50.  x2=50.  yl=.5  y2=3.675 
xcc=0.  ycc=0.  angla=60 
package  air 

rectangle  fill 

stations 

xs=0.0  yl=-8 . ,  0. 

end 
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A-2.  INPUT  DECK  FOR  THE  SECOND  CALCULATION 

A  typical  input  file  to  initiate  the  HULL  nm.  called  keel.in,  follows: 


keel  prob  00221.0200 
dimen=2  geom-1 
imax-50  jmax-100 
fail=l  strain=l 
nm=3  air=l  hha=2  walloy=3 
visc-1  dvisc-2 
nop“2  nstn=2 
header 

walloy  -  hha  0  1.1  km/s,  stochastic  stress+strain  failure 

mesh 

x0=-5  xmax=5  y0=-10  ymax=10 
generate 

package  walloy  v=l . leS 

rectangle  xl=-.414  x2=.414  yl=-10.26  y2=-.414 
circle  xc=0 .  yc=-0.414  radius^O . 414 
package  walloy  v=l . leS 

circle  xc=0 .  yc=-0.414  radius=0.414 
package  hha 

rectangle  xl=-50.  x2=50 .  yl=.5  y2=3.675 
xcc=0 .  ycc=0.  angla=60 
package  air 

rectangle  fill 

stations 

xs=0.0  yl—10.26,  0. 

end 
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